
############ linear system ###########

# declare 2 dependent variable ode system
# dx/dt = ax + by
# dy/dt = cx + dy
# input: 
#	t: current time point
#	var: c(x,y) current state
#	p: params, c(a,b,c,d)
# output: list(c(dx/dt,dy/dt))
linear.ode.system <- function(t, var, p) {
	x <- var[1]; y <- var[2]
	a <- p[1]; b <- p[2]; c <- p[3]; d <- p[4]

	return(list(c(
		a*x + b*y,
		c*x + d*y)))
}

# give x, find y such that 
# dx/dt = dy/dt = 0
linear.nullcline <- function(x, p) {
	a <- p[1]; b <- p[2]; c <- p[3]; d <- p[4]

	return(c(
		-a/b*x,
		-c/d*x))
}


########### run experiment ##########

	test.ode.system	= linear.ode.system 
	test.nullcline	= linear.nullcline 
	test.params		= c(-1,-1,-4,-5)
	test.xrange		= c(-2,2)
	test.yrange		= c(-2,2) 
	test.xstep		= 0.2
	test.ystep		= 0.2
	test.xlabel		= "x"
	test.ylabel		= "y"
	test.x.ini		= c(0.5)	
	test.y.ini		= c(-1.2)
	test.times		= seq(1,5,by=0.01)
	test.ode.root	= NULL
	test.ode.event	= NULL

draw.phase.plane(
			ode.sys = test.ode.system, ode.event = test.ode.event, ode.root = test.ode.root,
			nullcline = test.nullcline, 
			params = test.params,
			xrange = test.xrange, yrange = test.yrange, 
			xstep = test.xstep, ystep = test.ystep,
			xlabel = test.xlabel, ylabel = test.ylabel,
			x.ini = test.x.ini, y.ini = test.y.ini, times = test.times)

#draw.trajectory.over.time(
#						ode.sys = test.ode.system, ode.event = test.ode.event, ode.root = test.ode.root, 
#						params = test.params,
#						x.ini = test.x.ini, y.ini = test.y.ini, 
#						times = test.times,
#						xlabel = test.xlabel, ylabel = test.ylabel)

